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Abstract 

Most of numerical methods for deterministic simulations of rarefied gas flows use 
the discrete velocity (or discrete ordinate) approximation. In this approach, the kinetic 
equation is approximated with a global velocity grid. The grid must be large and fine 
enough to capture all the distribution functions, which is very expensive for high speed 
flows (like in hypersonic aerodynamics). In this article, we propose to use instead 
different velocity grids that are local in time and space: these grids dynamically adapt 
^.^ to the width of the distribution functions. The advantages and drawbacks of the 

\^ ', method are illustrated in several ID test cases. 
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1 Introduction 



X. 

^ . Most of deterministic numerical methods for rarefied flow simulations are based on a discrete 

- - -' velocity approximation of the Boltzmann equation, see for instance [20], [HI El El HI [151 [161 

In almost all these methods, the distribution function is approximated with a global 
velocity grid, for every point in the position space, for every time. This makes the method 
robust (conservation, entropy dissipation, positivity, stability, etc.) and relatively simple, 
but very expensive for many cases. Indeed, the grid must be large enough to contain all the 
distribution functions of the flow, and fine enough to capture every narrow distributions. 
The first constraint makes the grid very large for high speed flow with large temperatures. 
The second constraint makes the grid step very small, and hence a very large number of 
discrete velocities. This is for instance the case for atmospheric re-entry problems, where 
the flow is hypersonic. Such cases, especially in 3D, are very difficult to be simulated with 



such methods, due to the a discrete velocity grid that contains a prohibitively large number 
of points. 

Of course, particle solvers like the popular DSMC method do not suffer of such prob- 
lems [3]. However, if one is interested in deterministic Eulerian simulations, it is important 
to find a way to avoid the use of a too large number of discrete velocities. Up to our knowl- 
edge, there are a few papers on this subject. Aristov proposed in [1] an adaptive velocity 
grid for the ID shock structure calculation. However, the approach is very specific to this 
test case and has never been extended. More recently, Filbet and Rey |9] proposed to use a 
rescaling of the velocity variable to make the support of the distribution independent of the 
temperature and of the macroscopic velocity. Then the Boltzmann equation is transformed 
into a different form (with inertia terms due to the change of referential). In [2], Baranger 
et al. proposed an algorithm to locally refine the velocity grid wherever it is necessary and 
to coarsen it elsewhere. But this approach, which has been proved to be very efficient for 
steady flows, is still based on a global grid, and cannot be efficient enough for unsteady flows. 
Finally, Chen et al. [7] proposed to use a different velocity grid for each point in the position 
space and every time: from one point to another one, the grid is refined or coarsen by using 
an AMR technique. This seems to be very efficient, but all the grids have the same bounds 
(they all use the same background grid). It seems that a similar method was proposed at 
the same time by Kolobov et al., see [Ml [12]. 

In this paper, we propose a method that has several common features with the method 
of [22] and [9] but is still very different. The main difference is that each distribution is 
discretized on its one velocity grid: each grid has its own bounds and step that are evolved 
in time by using the conservation laws. The interaction between two space cells requires to 
use reconstruction techniques to approximate a distribution on different velocity grids. This 
paper is a preliminar work, for ID flows, that proposes a complete algorithmic approach. 
Several test cases illustrate the properties of the method and show its efficiency. 

The outline of the paper is the following. In section [2] is presented a simple ID kinetic 
BGK model and its standard discrete velocity approximation. In section |3l the local discrete 
velocity grid approach is described. The numerical tests are given in section HI 

2 A simple ID kinetic model and its standard velocity 
discretization 

We consider a one-dimensional gas described by the mass density of particles f{t,x,v) that 
at time t have the position x and the velocity v. The corresponding macroscopic quantities 
can be obtained by the moment vector U{t, x) = {mf{t, x, .)), where m(w) = (1, w, ^Iwp) and 
(0) = Jjj (j){v) dv for any velocity dependent function. This vector can be written component 
wise hy U = {p,pu,E), where p, pu, and E are the mass, momentum, and energy densities. 
The temperature T of the gas is defined by relation E = ^plu]"^ + \pRT, where R is the gas 
constant. 



The evolution of the gas is governed by the following BGK equation 

dJ + vdJ = -iMiU)-f), (1) 

r 

where M{U) is the local Maxwellian distribution defined through the macroscopic quantities 
f/ of / by 

and r = CT^ / p is the relaxation time. The constant C and u will be given in section |4] for 
each test case. 

From this equation, it is easy to establish the so called conservation laws that describe 
the time evolution of the moment vector U: 

dtU + d^ (mf) = 0. (3) 

For the numerical approximation of equation ([1]), a popular method is the discrete 
ordinate-or discrete velocity-method. It consists in choosing a grid V oi K points v^, and 
then in replacing the kinetic equation ([T]) by the finite set of K equations 

dtfk + vAh = -{Mk{U) - h) (4) 

r 

where /^(t, x) is an approximation of /(t, x, Vk) and MkiJJ) is an approximation of M{U){yk)- 
In order to describe the solution correctly, the discrete velocity grid V must capture 
all the distribution functions, that is to say at any time t, and for every position x. This 
means that V must be large enough to capture distributions with large mean velocity or 
large temperature, and fine enough to capture distributions with small temperature. See 
an illustration of this problem in figure [U In this figure, we show a 2D aerodynamical flow 
with three typical distributions functions (one in the upstream flow, one in the shock, and 
another one at the boundary). The corresponding velocity grid is shown in the same figure. 
Such a grid can be constructed as follows. We assume that a given distribution function 
whose corresponding mean velocity and temperature u and T can be well described if we 
restrict to the interval [u—AyBT^ u+AyBT]. This choice is reasonable, since in most physical 
problems, distributions are quite localized and they are concentrated in an interval which 
is not too far from the support of their local Maxwellian. Consequently, a first constraint, 
which necessary for all the distributions can be captured by the global velocity grid, is that 
its bounds Vmin and Vmax satisfy the following inequalities: 

Wmm < min(u(t,x) -4 v^MX^T^), Vmax > max{u{t,x) +4:y/RT{t,x)). (5) 



Since it is reasonable to require that there are at least three points between the inflexion 
points of any Maxwellian, the step of the global grid should be such that 



Av < min y/RT{t,x). (6) 

t,x 



Of course, such an approach requires to first estimate some bounds on the macroscopic fields 
that are global in time and space. Even if it is rarely mentioned in the literature, we believe 
that this is more or less the way people do when they use numerical simulations by using 
discrete velocity approximations. 

Note that the points of V are not necessarily uniformly reparted, since the grid could 
be refined only wherever it is necessary and coarsened elsewhere, as it is proposed in [8] for 
steady flows (with a simple and automatic way to define such a grid). However, for unsteady 
problems, the situation is more complex. Indeed, first, the estimations of the correct bounds 
and step of the grid are not necessarily available for every problems (the velocity or the 
temperature could reach much larger values that were not expected at some times of the 
simulation), like in complex shock interactions, for instance. Moreover, even if it is possible, 
this could lead to a grid which is extremely large and dense, hence leading to a very expensive 
simulation (see an example in section |1]2]). 

It is therefore very attractive to try to use local velocity discretization of the distribution 
function, which means to use a local discrete velocity grid (LDV) for each time and position. 
In other words, we would like to approximate f{t, x, .) with one different grid for each t and 
X. The clear advantage of this idea is that we can defined an optimally small grid for each 
distribution, thus we avoid the problems mentioned above. This approach is developed in 
the next section. 

3 A local discrete velocity grid approach 

In this section, we describe the basic algorithm that allows us to use local discrete velocity 
grids. Then we give some details about the reconstruction step, and we suggest possible 
definition of non symmetric grids. 

3.1 The method 

Now we assume that at time t„, the distribution function in each space cell [Xj_i,Xj_,_i] is 
approximated on a set V" of K local discrete velocities (the local discrete velocity grid, LDV 
for short) for i e {0,imax + !}• The cells [Xj_i,a;j_,_i] for i G {l;imax} represent the cells 
inside the domain. For simplicity, we assume here that all the local grids have the same 
number of points K, and the points are equally reparted. The first point is denoted by v^^i^ j 
and the last one by the w^^^, j. That is to say, we have for i G {0, imax + 1}: 

V? = H,k = <^n,^ + kAv^. k = 1 : K} , where A< 



K 



The velocity grids V^ and V^j^^x+i represent the velocity grids associated to the fictions cells. 
On this local grid, f(tn,Xi,.) is approximated by K values that are stored in the vector 
fr — i.fi'k)k=i- Each discrete value /";. is an approximation of f(tn,Xi,v^jJ. 



The corresponding macroscopic quantities U(tn,Xi) are approximated by Ul" with the 
following quadrature formula 



K 



U: = {mfnvr^ = E ^i<k)f",k^k, (7) 



fc=i 



where the Uk are the weights of the quadrature. 

When one wants to compute an approximation of / at the next time step tn+i, two 
problems occur. First, how to determine the local discrete velocity grid Vf"*"^? We will show 
below that this can be simply made by using the conservation laws. Second, how to exchange 
information between two space cells, since the local grids are not the same? This is where 
we use some interpolation procedure in the method. Let us now describe our algorithm step 
by step. 

Step 1: Macroscopic quantities at tn+i- 

We approximate the conservation relation ([3]) with the following first order upwind scheme 



* Ax \ "' 2 "2 



where the numerical fluxes are defined by 

$r+i = i^^^fDy. + (^^m/r+i>vn,, , (9) 



'i+i 



which is indeed an approximation of the flux (vmf{tn-,x^j^i)) at the cell interface. Here, 

we use the standard notation u ^ = (t; ± |t'|)/2. Note that each half flux is computed on the 
local velocity grid of the corresponding distribution. The vector [/" '* is an approximation 
of [/(t„+i, Xj), and we note u^^ '* and T""*" '* the corresponding velocity and temperature. 



n+l 



Step 2: discrete velocity grid V, 

We define this grid by using the new velocity and temperature ul '* and T" '* to get 
the bounds 



Then the new grid Vf """^ is defined as in the previous time step, that is to say by 

Vr' = {<r = <^n, + kAvr\ k = l:K}, with A<+^ = (v^;':^^^ - v-Jl,)/K. 

Step 3: distribution function at time t„+i. 

Here, equation ([1]) is approximated by a first order upwind scheme, with an implicit 
relaxation term. If the velocity variable is not discretized, we get for every v. 



At 



+ ^iM{ur'''){v)-f--^\v)). 



If now we take into account that each distribution /"^^, /", f^_i, and fl^_^i are defined on their 
own local velocity grid, this scheme must be modified by using a reconstruction procedure. 
Then we define the discrete values of /f^^ on its grid Vf ^^ by 



At 
7l' 



+ ^^,{M,{ur') - m') 



(11) 



for /c = 1 to K. In this relation, the distributions /", f^_ii and //^^ are used to reconstruct 
piecewise continuous in velocity functions ^"■, ^'1^, and f^j^^ that are defined as follows: 

■ <v <v'^ ■ 
fi{v) = { ^ , (12) 




where p" is a piecewise continuous function of v constructed through the values (Wj\, /f^)^;^, 
like a piecewise interpolated polynomial. This reconstruction will be discussed in next sec- 
tion. 

Step 4- Corrected macroscopic quantities. 

Here, we compute the moments of /"'''^ by the quadrature formula ([7]) to get 

K 

k=l 

Note that U^^'^ is an approximation of U{tn,Xi), like the vector [/" '* computed at step 1. 
This new step might be avoided by setting UJ^~^^ = f/f"*" '*, but it turns that it is required to 
ensure the positivity of our scheme, as it is stated below. 

Property 3.1. Assume that for every cell i, /" is non negative at each point of its local 
velocity grid and that the corresponding reconstructed piecewise function /" is non negative 
for every V. Finally, assume thatT^^ '* is positive for every i. Then under the CFL condition 
At < Ax/maxj^fedw^""^^!), /"^"^ is also non negative at each point of its local velocity grid, 
for every space cell. 

Proof. As it is standard for the upwind scheme for convection problems, it is sufficient to 
note that (TTT]) can be written as a linear combination of ^"', fl'_^, fi+ii ^^"^ M{U'^^ '*). The 
CFL condition of the proposition ensures that the coefficients of this combination are non 
negative, which gives the result. D 

To conclude this section, we comment on the assumption on the sign of T" '*. Of 
course, this temperature must be positive, since it is used to define the local velocity grid 



Vf"*"^ at step 2. However, it seems hardly possible to prove that this property is true for the 
scheme presented above. Nevertheless, we want to show below that a simple modification of 
the scheme is sufficient to ensure the positivity of this temperature. Indeed, we propose to 
replace the quadratures used to compute the macroscopic vector f/" and the numerical flux 
$" 1 (see (ITj) and ([9])) by the exact integral of the corresponding reconstructed functions, 

that is to say: 

If the reconstruction procedure uses a polynomial interpolation, these integrals are just 
half integrals of piecewise polynomial functions and can be evaluated explicitly. With this 
definition, we rewrite the discrete conservation law ([8]) as 

where </> is a piecewise continuous function of v. Now we have the following property. 

Property 3.2. Under the CFL condition At < Aa;/maxj.fc(|w"^.|), the function (p is non 
negative, and hence T^ '* is positive. 

Proof. Observe that is a linear combination of /f , fl^i, fi+i- The last two coefficients are 
always non negative. As a consequence of the CFL condition, the first coefficient, which is 
1 — ;^|f I, is non negative if \v\ is small enough (that is to say \v\ < maxj^fc(|w";i,|)). If v is 
larger, the coefficient is negative, but by construction f^iy) = 0. Consequently, is non 
negative for every v, and hence U^~^ '* is realized by a non negative distribution. It is then 
a standard result that the corresponding density, energy and temperature are positive. D 

Remark 3.1. Note that this CFL condition for the positivity of T" '* makes At depend 
on the v")^, while the CFL condition for the non-negativity of fl^~^^ (see property 13. ip makes 
At depend on the ff^^. This means that our modified scheme requires two time steps: a 
first one to compute f/" '*, and another one to compute fj^~^^. These time steps might be 
slightly different. 

However, we observe that in practice, even the original method (with quadratures and 
only one time step) preserves the positivity. This is why we do not use the modified scheme 
in the numerical tests presented in this paper. 

3.2 Interpolation step 

To compute fi{v'^^^) in equation fITT]) . we have to use a reconstruction procedure. First, 
if v^~^^ is external to Vf , we set /^(f"^^) = 0. If w"^^ is inside Vf , it is not a node of 
Vf in general, and we use a piecewise polynomial interpolation. We observed that first 
order polynomial interpolation is not accurate enough. However, higher order interpolation 



produces oscillations, especially in very rarefied regimes, which is probably due to the large 
velocity gradients (an even discontinuities) of the distribution functions in such regimes. 

Consequently, we use the essentially non oscillatory (ENO) interpolation (see [19] or [TO]): 
with 3 or 4 point ENO interpolation, the accuracy is good and there is almost no oscillation. 

The reconstruction algorithm is summarized below. 

1. If K+^l > max(te,„,|, te,,,|), then set /;"«+^) = 

2. else 

(a) find the interval K"fc/,'?^rfc/+i] inside Vf that contains w"^"^ 

(b) compute the q point ENO polynomial function P defined on the stencil {f f;./_( ^^n, v'iku^q} 



(c) set frKt') = Pi<^'' 



Of course, the same procedure is applied to the other reconstructed values /j'!^i(ff^^) and 
flUiv-^') in dUD. 

3.3 Non symmetric local discrete velocity grids 

When the flow is far from equilibrium, the distribution function are different from their local 
Maxwellian, and might have a non symmetric shape. In particular, their support might be 
non symmetric as well (see for instance the heat transfer problem in the rarefied regime as 
shown in section 14.31 However, the local grids defined in section 13.11 are based on the local 
Maxwellian are are necessarily symmetric. In this section, we propose a method to modify 
the grid if necessary. 

First, note that up to now, we have defined uniform local grids with a constant number 
of points. However, the method is readily extended to non uniform grids with a variable 
number of points : we just have to replace K by K^ ad Uk by a;";, in every expressions of 
section 13.11 

Then, the idea is to enlarge the grid Vf"*"^ if /"^^ is not small enough at its boundaries. 
This is made as follows. 

We perform a splitting between the relaxation step and the transport step. We first 
compute the intermediate quantity /"^ by using the transport equation as 



(14) 



for k = 1 to K. At the end of the transport step, the values of the distribution function 
/f at the boundary points ffmax ^^d v^^^ of Vf"*"^ are compared to the maximum value 
of the distribution on the grid. If these boundary values are not small enough, new grid 
points wr = ff^Q^ + At^J^"*"^ and wl = v'^'^m ~ ^'^^l^^ are added outside the grid by using 
again the transport step, until they are small enough. At the end of this step, the modified 



grid Vf"*"^ now has K^~^^ velocities. Finally, /"^^ is obtained from the relaxation step through 
the relation 

fS' = fT' + -^(M.iUr') - f-n (15) 

i 

for A; = 1 to Kf+^ 

It would also be interesting to use an automatic refinement of the local grid around the 
possible discontinuities, but this is not studied in this paper. 

4 Numerical results 

In this section we present three numerical tests to illustrate the main features of our method 
(denoted by LDV, for local discrete velocity grid). It is compared to a standard discrete 
velocity method (with a global velocity grid) denoted by DVM. First, the numerical scheme 
is tested on the Sod test case for three different regimes: the rarefied, fluid, and free transport 
regimes. The second test is the two interacting blast waves problem in which very high 
temperature differences make the standard DVM unefficient. The third test case is devoted 
to the heat transfer problem, where the use of non symmetric local grids is shown to be 
necessary when the Knudsen number increases. In these tests, the gas constant R is 208.1, 
except for the free transport regime in section 14.1.31 where R= 1. 

4.1 Sod test case 

4.1.1 Rarefied regime 

We consider a classical Sod test in a rarefied regime for the BGK model ([1]) with the param- 
eters uj = —0.19 and C = 1.08 ■ 10"^ used in the relaxation time r. The space domain is 
the interval [0, 0.6] discretized with 300 points. The initial condition is the local Maxwellian 
distribution whose macroscopic quantities are given by: 

T(x) = 0.00480208, p(x) = 0.0001, u{x) = for x G [0, 0.3] 

T(x) = 0.00384167, p{x) = 0.0000125, u{x) = for a; g]0.3, 0.6]. ^ '* 

The DVM and LDV methods are compared to a reference solution given by the DVM 
method computed for a large and fine velocity grid (obtained after a convergence study). 
This velocity grid has 600 points uniformly reparted in the interval [vmin,v„iax] where 

Vmin =mm{u(t,x) - 6^/RT{t,x)), Vmax = max{u{t, x) + 6^/RT{t, x)) (17) 

Of course, such bounds cannot be determined a priori, and several computations with larger 
and larger velocity grids have to be done before the correct bounds are found. This illustrates 
the difficulty to use a standard DVM when the extreme values of the velocity and the 

9 



temperature are not known a priori. Indeed, here, the temperature in the shock after the 
initial time is higher than the two initial left and right temperatures. If the velocity grid is 
computed with formula ([5]) and ([6]), when the bounds are estimated with the initial values 
of T and u, then it is not large enough, and the results are not correct, even if the number 
of velocities is increased so as to reach the grid convergence which is obtained here with 100 
velocities. This is shown in figure |2l 

At the contrary, our LDV method dynamically adapts to the time variations of u and 
T and gives very accurate results with 30 discrete velocities only in each local grid, as it is 
shown in figure O Consequently, the LDV method is very efficient for this case. 

Note that in this test, the reconstruction procedure used the 4-points ENO interpolation, 
see section 14.1.31 for an analysis of the influence of the order of interpolation. 

4.1.2 Fluid regime 

Now, we consider the same Sod test case, but in the fluid regime. This regime corresponds 
to the limit case of equation ([1]) when r is set to 0. Note that since both DVM and LDV 
methods are used with a time semi-implicit scheme, taking r = means that Z""*"^ is set 
to M{U'^~^^) at each time step, and hence we get two different numerical schemes for the 
compressible Euler equations of gas dynamics. 

Here, the reference DVM has 100 velocities only. The DVM with the grid computed with 
the initial temperatures and velocity has also 100 velocities, but it still gives incorrect values 
(see flgure[3]), for the same reasons as mentioned in the previous section. At the contrary, 
our LDV method gives very accurate results with 10 velocities only. 

Note that here, step 3 of our scheme is nothing but the projection step fl^~^^ = M{u^~^ '*), 
and hence the choice of the interpolation procedure has no influence. 

4.1.3 Free transport regime 

We consider the free transport regime corresponding to the limit case in ([1]) when r tends to 
+00. In this subsection, we take R = 1 and the standard non dimension values for the Sod 
test case. The space domain is [—1, 1] and is discretized with 300 points. The distribution 
function is initialized by the local Maxwellian distribution function with 

p= 1, M = 0, p = 1, in [-1,0[, 
p = 0.125, n = 0, p = 0.1in [0,1]. ^ ^ 

For this test case the numerical results are compared to an analytical solution of the free 
transport equation. 

It is very difficult to accurately approximate the free transport equation with a stan- 
dard DVM: the macroscopic profiles obtained with the DVM show several plateaux. These 
plateaux are not due to the space and time approximation, but are only due to the velocity 
discretization. Indeed, it can be easily proved that the macroscopic profiles of the DVM so- 
lution at time t have plateaux of length tAv, where Av is the step of the global uniform grid. 

10 



This is clearly seen in figures |1H6] (top), where the DVM has 30 discrete velocities. This phe- 
nomenon is known as the ray effect that appears with the discrete ordinate approximations 
of the radiative transfer equation. 

When we test this problem with our LDV method, with 2 point piecewise interpolation, 
and 30 velocities in each local grids, the results are very bad: we observe very large oscillations 
(see figures HHHl bottom). This is probably due to the fact that this interpolation is not 
accurate enough. Then we use 3 and 4 point ENO interpolation, and we observe in the same 
figures that the solution is now much closer to the exact solution. Moreover, while we have 
the same number of discrete velocities in each local grid as in the DVM, we observe that the 
plateaux are completely eliminated. 

However, there are still some oscillations in the results obtained with the LDV. This is 
probably due to the fact that these results were obtained without the moment correction 
step 4 of our algorithm (see (fT3|) ). Indeed, if we do now the same simulation with this 
moment correction step, the results are very good: there are no oscillation at all, no plateau 
phenomenon, and the results are much closer to the analytical solution, see figures U\M 

4.2 Two interacting blast waves 

This section is devoted to the test case called "the two interacting blast waves" (see [T7]). 
Here, the relaxation time is defined with u = —0.19 and C = 1.08 10"^. The space domain is 
the interval [0, 1] which is still discretized with 300 points. The initial distribution function 
is a local Maxwellian distribution which macroscopic quantities are given by p = 1 and u = 
everywhere, and 

T = 4.8, in [0, 0.1], T = 4.8 10"^ in ]0.1, 0.9], T = 4.8 10"^ in ]0.9, 1]. 

On the left and right boundaries, we use Neumann boundary conditions. 

Here, the bounds of the global grid of the DVM are determined by the largest initial 
temperature, and its step size is given by the smallest initial temperature. Then we find 
that the coarsest global grid that satisfies conditions ([3]-E]) has not less than 2 551 velocities! 
In figures [T0l - [T2t we observe that the LDV method requires only 30 velocities to give results 
that are very close to the DVM method (with 2551 velocities), both before and after the 
choc. This proves the high efficiency of the LDV approach for this case. Note that here 
again, we use a 4 point ENO interpolation in the LDV method. 

4.3 Heat transfer problem 

In this test, we consider the evolution of a gas enclosed between two walls kept at temperature 
Tl = 300 and Tr = 1000. At these walls, the distribution function satisfies the diffuse 
boundary condition 

/(a; = 0,i;>0) = M,,,o,T„ /(x = 1, i; < 0) = M,^,o,t, (19) 
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where 



^^ _ L<oVf{^ = 0,v)dv ^^ ^ l^^vf{x = l,v)dv 



and Mpfl^T denotes p/v27rRTexp{—v'^/2RT) for every p and T. The initial data is the 
Maxwelhan with density po (to be defined later), velocity uq = 0, and Tq = Tl = 300. Here, 
the relaxation time is defined with u = —0.5 and c = 6.15 ■ 10~^. 

The boundary conditions are taken into account in our numerical scheme by a ghost cell 
technique, as it is standard in finite volume schemes. Left and right ghost cells are defined for 
i = and i = imax + 1, and the velocity grids V^ and Vj*^^^.,.]^ in these cells are defined so as 
to correctly describe the corresponding wall Maxwellians. Then, the density pl and p/j, that 
are defined as the ratio of an outgoing mass fiuxes at a wall to the corresponding incoming 
Maxwellian mass fiux, are approximated by using the boundary cell and the corresponding 
ghost cell, that is to say: 

^^" (^;+M(l,0,r^))vy' ^''~ {v-Mil,0,Tn))vr ^, ' 

For this test case several Knudsen number can be used. This number is parametrized by 
the initial density po- We first analyze the LDV method in the transitional regime (po = 0.1, 
which gives Kn = 10~^): here, both the LDV and DVM are converged with 30 velocities, 
while the space domain [0, 1] is discretized with 1000 points. However, the results obtained 
with the LDV are not accurate enough (see figure [T3l) . An analysis of this problem shows 
that this is due to the local grid close to the right boundary which is not large enough: for 
small times, the distributions at these points are highly non symmetric. 

To correct this problem, we use the algorithm proposed in section [3731 to enlarge the local 
grids in a non symmetric way. Then, the LDV with 30 velocities now gives results that are 
indistinguishable from the DVM, see figure [TH 

Finally, we test the LDV method in the rarefied regime {Kn = 1): the LDV (with 
enlarged non symmetric local grids) and the DVM are converged with 300 velocities, while 
the space domain [0, 1] is discretized with 300 points. Here again, both methods give results 
that are almost indistinguishable (see figure [13]) . Unfortunately, the number of velocities 
required to get converged results is very large here, even for the LDV (300 velocities). This 
is probably due to the fact that the distribution function is discontinuous with respect to 
the velocity in this test, with very large jumps: our velocity grids, even the local ones, are 
uniform, and cannot capture these discontinuities when the number of velocities is too small. 
However, we show in figure dn] the results obtained with 50 velocities only, and we observe 
that the LDV gives results that are more accurate than the DVM. 

Remark 4.1. Note that in these tests on the heat transfer problem, we use the moment 
correction step, otherwise, some sharp oscillations appear at some points. For large Knudsen 
numbers, this correction step is essential, since the oscillations can be very large. 

Moreover, note the interpolation procedure has been performed with a 4 points ENO 
method, and no improvements are observed if we use larger stencils. 
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5 Conclusion and perspectives 

We have presented a new velocity discretization of kinetic equations of rarefied gas dynamics: 
in this method, the distribution functions are discretized with velocity grids that are local 
in time and space, contrary to standard discrete velocity or discrete ordinates methods. The 
local grids dynamically adapt to time and space variations of the support of the distribution 
function, by using the conservation laws. 

This method is very efficient in case of strong variations of the temperature, for which a 
standard discrete velocity method requires a very large number of velocities. Moreover, it 
allows to eliminate the plateau phenomenon in very rarefied regimes. 

The extension of this method to 2D problems is in progress and will be presented in a 
forthcoming work. 

However, several aspects of the method should be better understood, in particular, why 
are there some oscillations if the moment correction step is not used, in the rarefied and free 
transport regimes, even with high order ENO interpolation? A mathematical analysis of the 
numerical method could be interesting here. 
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Figure 1: Three distribution functions in different space points of a computational domain 
for a re-entry problem, and the corresponding global discrete velocity grid V. 
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Figure 2: Sod test case, rarefied regime: temperature (top) and velocity (bottom) profiles, 
at time 7.34 10~^. The solid line is reference solution obtained with the DVM with the 
enlarged global grid (600 points), the dotted line is the DVM with an incorrect grid (100 
points), while the dot-dashed line is the LDV method (30 points). 
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Figure 3: Sod test case, fluid regime: temperature (top) and velocity (bottom) profiles, at 
time 7.34 10~^. The solid line is reference solution obtained with the DVM with the enlarged 
global grid (100 points), the dotted line is the DVM with an incorrect grid (100 points), while 
the dot-dashed line is the LDV method (10 points). 
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Figure 4: Sod test case, free transport: temperature at time 0.3. Top: comparison between 
the exact solution (dot-dashed) and with a 30 points DVM (sohd). Bottom, comparison 
between the exact solution (dot-dashed) and several LDVs with 30 points: without ENO 
('o'), with 3 points-ENO (dotted), with 4 points-ENO (solid). 
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Figure 5: Same as figure H] for the velocity. 
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Figure 6: Same as figure H] for tlie density. 
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Figure 7: Sod test case, free transport: temperature at time 0.3, comparison between the 
exact solution (dot-dashed) and several LDVs with 30 points and the moment correction 
step: without ENO ('o'), with 3 points-ENO (dotted), with 4 points-ENO (solid). 
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Figure 8: Sod test case, free transport: velocity at time 0.3, comparison between the exact 
solution (dot-dashed) and several LDVs with 30 points and the moment correction step: 
without ENO ('o'), with 3 points-ENO (dotted), with 4 points-ENO (sohd). 
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Figure 9: Sod test case, free transport: density at time 0.3, comparison between the exact 
solution (dot-dashed) and several LDVs with 30 points and the moment correction step: 
without ENO ('o'), with 3 points-ENO (dotted), with 4 points-ENO (sohd). 
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Figure 10: "Two interacting blast waves": Temperature, before the choc at time 0.008 (top) 
and after the choc at time 0.02 (bottom). The sohd hne is the solution obtained with the 
LDV method (30 points), the dotted line is the solution computed with the DVM method 
(2 551 velocities). 
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Figure 11: "Two interacting blast waves": velocity (same as figure [TO 
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Figure 12: "Two interacting blast waves": pressure (same as figure [TO 
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Figure 13: Heat transfer problem, transitional regime: temperature (top) and velocity (bot- 
tom) at time 1.3 10~^, Kn = 10~^. The solid line is the solution given by the LDV method 
(30 velocities), the dotted line is the DVM method (30 velocities). 



27 



900 



800 



700 



600 - 



500 - 



400 



300 





-50 L 



0,1 



0,2 



0,3 



0,4 



0,5 



0,6 



Figure 14: Heat transfer problem, transitional regime (test of the non symmetric local grid): 
temperature (top) and velocity (bottom) at time 1.3 10~^, Kn = 10~^. The solid line is the 
solution given by the LDV method (30 velocities), the dotted line is the DVM method (30 
velocities). 
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Figure 15: Heat transfer problem, rarefied regime: temperature (top) and velocity (bottom) 
at time 1.3 10~^, Kn = 1. The solid line is the solution given by the LDV method (300 
velocities, non symmetric local grids), the dotted line is the DVM method (300 velocities). 
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Figure 16: Heat transfer problem, rarefied regime: temperature (top) and velocity (bottom) 
at time 1.310^'^, Kn = 1. The solid line is the solution given by the LDV method (50 
velocities, non symmetric local grids), the dotted line is the DVM method (50 velocities), 
the dashed line is the DVM method with 300 velocities (reference solution). 
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